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ABSTRACT 


Nunerical  integration  methods  for  solution  of  the  system  of 
differential  equations  found  in  ballistic  rocket  trajectory  programs 
are  discussed. 

The  general  discussion  entails  the  explicit  formulas  of  Runge- 
Kutta  and  predictor-corrector  methods  and  their  errors,  and  a  brief 
description  of  other  methods  that  could  be  employed. 
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INTRODUCTION 


In  a  ballistic  rocket  trajectory  simulation  program  the  system 
of  differential  equations  used  to  describe  the  ballistic  model  is 
a  highly  complex  system.  In  particular  the  six -degree  of  freedom 
model  used  by  the  Atmospheric  Sciences  Office  at  White  Sands  Missile 
Range,  consists  of  a  system  of  twenty-one  2nd  order  ordinary  differen¬ 
tial  equations  which  are  to  be  solved  for  the  ballistic  rocket's  com¬ 
ponents  of  acceleration,  velocity,  and  position  at  discrete  time  in¬ 
tervals.  The  most  feasible  method  for  solving  the  system  is  to  pro¬ 
gram  it  for  a  computer  and  solve  by  some  numerical  integration 
procedure. 

This  report  will  present  several  numerical  integration  schemes 
which  are  currently  being  used  or  are  feasible  for  use  in  a  ballistic 
rocket  simulation  program.* 


RUNGE-KUTTA  METHODS 
Formulas 


Consider  the  system  of  first  order  differential  equations 
d/i 

yJ  “  33^-  =  yj(x),  y2(x), 


(1) 


satisfying  the  initial  conditions 


y.  (x  )  =  y.  . 
o'  'lO 


(2) 


We  desire  the  values  of  y^^Cx^  ♦  h),  i  »  1,  2,  ...,  N,  where  h  is  the 
increment  in  the  independent  variable. 


•This  report  was  presented  to  the  Second  Annual  Conference  on  Pure 
and  Applied  Mathematics  held  at  New  Mexico  Institute  of  Mining  and 
Technology,  Socorro  New  Mexico,  feb.  24-25,  1967. 


The  Rungo-Kutta  method  is  an  algorithm  designed  to  approximate 
the  Ta/lor  series  solution 


y.  (x  +  h)  »  y.  ♦  h  y 


(1) 


(2) 


+  . . . 


(3) 


of  (1),  where  (*(,)  denotes  the  kth  derivative  and  i  =  1,  2, 

N.  However,  in  deriving  the  defining  equations  of  the  algorithm, 
the  use  of  the  appropriate  assumptions  makes  it  unnecessary  to  eval¬ 
uate  and  define  the  derivatives  of  higher  order. 


The  general  equations  defining  a  fourth  order  Runge-Kutta  method 
are  given  by 

=  hf(x^  ♦  ah,  y^  + 

•^3  *  *  “l^’  ^0  ^  ^1*^1 

=  hf(x^  +  02h,  y^  +  62^^  +  Y2*^2  *  ^2*^3^ 

K  =  y(x^  +  h)  -  y^  =  Pjkj  +  112*^2  *  ^3*'3  ^4*^4 


where  the  bar  denotes  an  approximate  value. 

The  thirteen  parameters  involved  are  found  by  expanding  K  about 

the  point  (x  ,  y  )  so  as  to  agree  with  the  Taylor  series  expansion 
00  ^ 

(3)  up  to  and  including  terms  of  order  h  ;  thus,  the  fourth  order 
method.  Initially,  the  solution  will  result  in  eight  equations  in 
thirteen  unknowns.  However,  under  the  appropriate  apsumptions,  the 
system  will  resolve  into  a  system  of  eight  equations  in  nine  vnknowns, 
hence  a  one  parameter  family. 
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One  solution  is  the  classical  formulas  of  Kunge,  given  by. 


*^1  *  ^^^*0’  ^0^ 

^2  -  hf(x^  -  T-  *  r’ 

Kj  =  hf(x^  •*  1.  Xo  *  '5) 

=  hf(x^  +  h,  y^  t  Kj) 

I  -  ^  CKi  *  2K2  ♦  21Cj  *  K^) 


One  will  notice  that  when  y'  =  f(x),  Runge's  formulas  reduce  to  Simp¬ 
son's  rule: 

X  ♦h  ,  , 

K  =  I  °  f(x)dx  =  I  [f(x^)  ♦  4f  (x^  ♦  |)  4  f(x^  4  h)]  (6) 

*0 


Another  solution  to  the  one  parameter  family  is  the  classical 
formulas  of  Kutta  given  by 


*^1  =  ^0^ 

k2  =  hf(x^  *  p  Xo  * 

2h  ’^l  *^2 

X3  =  hf(x^  *  r.  Xo  -  r  *  r  ^ 

»=  hf(x^  4  h,  Xj,  +  Kj  -  4  Kj) 

K  ■=  j  (Kj  ♦  3K2  4  SKj  4  K^) 


(7) 
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To  apply  the  Runge-Kutta  method  using  a  high  speed  digital  com¬ 
puter,  Gill  [8]  developed  a  calculation  procedure  which 

(a)  requires  a  minimum  number  of  storage  registers, 

(b)  gives  a  high  accuracy, 

(c)  requires  comparatively  few  instructions. 

Gill  accomplished  the  above  by  defining  an  auxiliary  set  of  pa¬ 
rameters  in  the  initial  solution  set,  thereby  obtaining  the  following 
solution; 

Kl  •  hfCx,,.  y^)  Xl  ■  y.  * 

h  ■  yp  y^  •  yj  ♦  (1  -  )  (Kj  -  Qj) 

Kj .  hf(x^ .  yp  y,  ■  ya  *  <1  *  )  ‘S  ■  '’2’ 

•^4  “  hf(x^  *  Xs)  y4  *  Xq  ^  “  y 

K  «  Ki  +  j  (1  -  /j  )  K2  ♦  j  (1  ♦  /  j  )  Kj  ♦  I  K4  (8) 

Qj  -  ♦  3[1/2(K^  -  2Q^)]  -  ^  Kj 

Q2  '  Qi  *  3[(1  -  /I  )  (K^  -  Qj)]  -  (1  -  )  ^2 

Q3  .  9-  ♦  3[(1  ♦  /I  )  (K3  -  Q2)1  -  (1  +  /|  )  K3 

Q4  »  Q3  *  3[i  (K4  -  2Q3)]  -  i  K4 
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To  obtain  a  general  expression  for  truncation  error  per  step, 
the  functions  1(  =  y(x  ♦  h)  -  y  and  K  =  y(x^  +h)  -  y^  are  expanded 

through  terms  of  h^,  thus  giving  the  per-step  truncation  error  as 

E  ^  r  - 

♦  3f  P  ♦  2f  ST  -  4f^  T)  ♦  2f  T^]  ♦  . . .  (9) 

yy  y  y  y 

where  T*^  =  o’^f,  =  o'^fy.P  =  (Df)^. 

Lotkin  [12]  determined  a  bound  for  L  which  is  given  by 

|E|  !  mlV  (10) 

where  in  R*,  a  region  containing  (x^,  y^), 

|f(x,  y)  I  !  M. 


(11) 


where  the  positive  constants  M  and  L,  are  independent  of  (x,  y), 
for  i  ♦  j  ^  4. 
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A  nuaerical  integration  method  is  called  stable  if  at  the  nth 
step,  the  total  error,  which  is  both  truncation  and  round-off  error, 
is  at  a  ninimuou  thereby  forcing  the  approximate  solution  to  tend  to 
the  true  solution. 


Hence,  as  for  the  stability  of  the  Runge-K^tta  methods,  Carr 
[3]  proved  a  theorem,  of  which  the  essential  parts  are  stated  below: 

Let  f^  be  continuous,  negative  and  bounded  from  above  and 
below  in  some  region  R  of  the  (x,  y)  plane,  i.e.,  >  0,  fy 


e(-M„ 


M^);  further,  let  E  be  the  absolute  value  of  the  maximum 


error  introduced  at  each  step,  and  0*  be  a  region  in  which  the  so¬ 
lution  of  the  difference  equation  tends  to  the  y-boundary  of  D  no 
closer  than  Qh  ^  | e. j ,  where 


K,  K, 

Q  >  max  |f(x,  y)  |  i  max  [|^|,|^|,  |k.1]  (12) 

x,yeD  2 


and  e.  is  the  error  at  the  ith  step.  Then,  the  Kutta  fourth  order 

method  has  a  bound  on  the  total  error  in  the  ith  step,  where  this 
bound  is  given  by 


(13) 


in  some  region  D*,  and  where  the  step-size,  h,  is  given  by 


h  <  min  ( 


(14) 


If,  -  0  (the  unstable  case),  but  bounded  throughout  U, 
5  fy  f  M,  then  the  propagated  error,  in  U*  at  the 

(i  +  1)  st  step  is  given  by 


X 


hM 

e 


(15) 
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and 


) 


and  E  and  h  are  as  given  before. 

Statements  (13),  (14),  (15)  relate  the  step-size  and  the  prop¬ 
agated  error.  Thus,  an  h  can  be  found  that  will  make  the  propagated 
error  less  than  a  certain  bound,  if  this  bound  is  the  known  bound 
on  the  partial  derivative.  Algorithms  for  finding  such  step«6izes 
can  be  found  in  Carr  [3] . 

Although  various  boiutds  on  the  truncation  error  are  known,  such 
as  the  one  given  above,  their  usage  in  a  computer  program  is  usually 
not  feasible.  Thus,  some  practical  scheme  for  estimation  of  this 
error  must  be  devised. 


One  such  scheme  was  devised  by  Richardson  [5] .  This  scheme 
is  based  on  numerically  integrating  with  step-sizes  of  h  and  k  and 


comparing  the  results  using  these  step-sizes.  A  variation  of  this 
which  we  use  is  the  following: 


.(1) 


the  value  ob- 


Let  Y  be  the  true  value  of  y  at  x  ♦h,  y 

(2)  ° 

tained  at  x  ♦h  using  hj  »  h,  y'^  the  value  obtained  at  x,+h 


using  h~  =  h  .  Then  for  small  h, 
*  0 


^  „(2) 


The  scheme  in  (16)  can  be  used  to  check  the  validity  of  the  answer 
for  the  purpose  of  halving  and  doubling  the  step-size  h. 


(16) 


Analysis 


By  virtue  of  the  Runge-Kutta  methods  requiring  only  information 
from  one  previous  step,  the  method  has  desirable  stability  charac¬ 
teristics  and  ease  of  halving  or  doubling  the  step-size  h. 

It  has  been  shown,  however,  by  Blum  [2]  and  Fyfe  [6],  that  through 
proper  redefinition  of  the  solution  set  of  parameters,  most  other 
versions  of  the  fourth  order  Runge-Kutta  method  can  avail  themselves 
of  the  reduced  storage  of  the  Gill  algorithm. 
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If  the  first  derivative  of  the  functions  are  very  involved,  then 
several  evaluations  of  these  first  derivatives  may  be  somewhat  time 
consuming  and  thus  cause  the  method  to  be  uneconomical  at  times. 
Another  distinct  shortcoming  of  the  method  is  that  of  the  error. 
Neither  truncation  error,  nor  its  estimate,  is  obtained  in  the  cal¬ 
culation  procedure,  thus  necessitating  the  approximate  comparisons 
described  previously. 

Finally,  the  Runge-Kutta  methods  can  be  used  as  a  starting  pro¬ 
cedure  for  other  methods,  such  as  multistep  methods. 


MULTISTEP  METHODS 

Standard  Predictor -Corrector  Methods 


A  multistep  method  is  a  method  in  which  the  calculation  of  the 
y^'s  at  the  (M+1)  st  step  depends  on  knowledge  of  the  yi's  and  the 

f^'s  at  the  Mth,  (M-1)  st,  (M-2)nd,  etc.  steps.  A  method  is  called 

aliTm-step  method  if  only  m-steps  of  previous  information  are  required. 
In  the  following  discussion  only  four-step  methods  will  be  considered 
since  methods  with  fewer  back  values  can  be  extracted  from  the  dis¬ 
cussion  below. 

The  multistep  method  with  which  we  are  concerned  here  is  known 
as  the  predictor-corrector  method.  This  method  requires  a  formula 
for  finding  a  first  estimate  of  each  y^^  (hence,  the  predictor);  and 

then  evaluating  each  function  f^  we  substitute  this  into  a  formula 

which  will  adjust  the  value  obtained  from  the  predictor,  hence  the 
corrector.  By  a  standard  method  we  mean  one  in  which  the  same  step- 
size  is  used  in  all  the  equations  in  integrating  for  each  y^^. 

First  of  all,  rewrite  the  system  of  differential  equations  in 
(1)  as 


^i'  “  ^i^^’  ^1’  ^2’  •••’ 


(17) 


with  the  initial  conditions 


yi(Xo)  -  yi(o) 


(18) 
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Ke  desire  the  values  y^(x^  ♦h),  i  >  1,  2,  ...  N  where  h  is  again  the 
increment  in  the  independent  variable. 

Next  let  us  clarify  some  notation  which  we  will  be  using. 

y^CM)  »  y^^  (x^  ♦Mh)  i  «  1,  2,  ....  N 

P.(M+1)  -  y.(x^  ♦(M+l)h)  i  »  1.  2,  ...,  N  (19) 

y^'CM)  •  yj’(xQ+Mh)  »  f^(x^+Mh,  y^CM),  y2(M),  ....  y^jCM)) 

7.'(M*1)  »  y.'(x^+(M+l)h)  =  f.(x^*CM*l)h.  Pj(M*l).  ...,  Pjj(M*l) 
C.(M*1)  =  y.(Xp^(M+l)h)  i  <=  1,  2,  ....  N 

By  using  (17).  the  defining  equations  for  the  standard  predictor- 
corrector  algorithm  are  given  by 

P.(M*1)  «  ajy.(M)  ♦  bjy.(M.l)  ♦  Cjy.(M-2) 

♦djy.(M-3)^h[ejy.'(M)  +  fiy.'(M-l)  (20) 

^  giyi'(M-2)  ♦  Kjy.'(M.3)l,  i  «  1.2,  ....  N 
C.(M^l)  =  a2y. (M)  ♦  b2y.(M.l)  ♦  C2y.(M-2) 

+  h[d2y.'(M+l)  ♦  e2y.'(M)  ♦  f2yi’(M-l)  (21) 

♦  g2yi'(M-2)l.  i  »1,  2,  ....  N 

where  (20)  is  the  predictor.  (21)  the  corrector,  and  the  coefficients 
are  to  be  determined. 
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Using  the  standard  technique  of  [10]  for  solution  of  the  con¬ 
stant  coefficients,  the  solution  set  is: 


Sj  ■  9  -  dj  -  3ej  +  3Kj 


bj  »  9  -  9dj  ♦  24Kj 


Cj  «  -  17  ♦  9dj  ♦  3ej  -  27Kj 


a2  “  9  -  15d2  -  3e2 


b2  »  9  -  24d2 


c_  =  -  17  ♦  39d»  +  3e, 


dj  =  dj 


d2  =  d2 


(22) 


®1  “  ®1 


®2  “  ®2 


fj  =  -  18  ♦  6dj  +  4ej  -  17Kj;  f2  *  -  18  ♦  39d2  ♦  4e2 
gj  =  -  6  +  6dj  ♦  Bj  -  14Kj;  *  -  6  ♦  14d2  ♦  62 


Ki  =  Ki 


Since  the  coefficients  in  (22)  are  in  terms  of  the  parameters 
^1’  ®1*  *^1’  ^2*  ®2*  constants  can  be  determined  so  as  to  give 

desirable  stability  characteristics. 

The  truncation  error  in  the  predictor-corrector  equations  is 
given  by 


Ep  =  (9  .  Jdj  -  ej  -  lOKj)  (23) 

=  30 


where  and  are  the  predictor  and  the  corrector  truncation  err^r. 

respectively.  However,  the  predictor  and  the  corrector  were  chosen 
independently  of  fourth  order.  Consequently,  from  Henrici  [11  ,  p 
261), the  truncation  error  of  the  algorithm  is  given  by  E^  alone. 
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But,  may  not  be  available;  thus,  s<»ie  means  for  reflecting  the 
truncation  error  must  be  made  available. 

The  standard  technique  for  adjusting  P^(M+1)  with  respect  to 
the  local  truncation  error  is  to  add  . 

E 

E(M+1)  =  (P.(M)  -  C.(M))  (24) 

C  p  V 


to  P,(M+1),  to  improve  its  value.  Notice  the  fact  that  E(M+1)  will 

^  fSl 

not  involve  y  . 

Similarly,  the  common  method  of  keeping  track  of  the  total  trun¬ 
cation  error  is  to  use 

E 

E^(M+l)  =  (P.(M+1)  -  C  (M+D)  (25) 

c  p 


as  an  estimate  of  this  error  at  the  (M'«-l)  st  step,  and  from  (25) 
criteria  for  halving  and  doubling  the  step^ize  h  can  be  generated. 

According  to  Crane  and  Klopfenstein  [4],  the  values  for  the  con¬ 
stants  in  (22)  which  yield  a  stability  situation  not  unlike  that  of 
Runge-Kutta  methods  are  given  by  ' '' 


a^  =  1.54765200 

bj  =  -1.86750300 
Cj  =  2.01720400 
d^  *  -0.697353000 
e^  »  2.00224700 
fj  9  -2.03169000 
gj  =  1.81860900 
Kj  *  -0.714320000 


3.2  •  1-0000 

b^  «  0. 

C2  *  0* 

d^  »  0.375000000 
e^  *  0.791666667 
£2  *-0.208333333 
g2  *  0.0416666667 


(26) 
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where  the  trailing  zeros  are  used  tc  make  the  method  as  "numerically" 
fourth  order  as  possible.  The  coefficients  in  the  corrector  equation 
determine  the  standard  Adams  -  Moulton  corrector. 


The  Generalized  Predictor-Corrector 


By  using  the  standard  predictor-corrector  formulas  as  discussed 
in  the  preceding  section,  one  inherent  shortcoming  arises;  that  of 
using  the  same  fixed  step-size  h  in  integration  of  the  various  y^^'s. 

In  many  physical  situations,  and  particularly  in  the  mathematical 
model  representing  the  trajectory  of  a  ballistic  rocket,  the  system 
of  differential  equations  that  arise^'may  have  some  solutions  vary¬ 
ing  much  more  rapidly  than  others  with  respect  to  the  independent 
variable.  Since  the  step-size  must  be  chosen  to  give  the  desired 
accuracy  in  the  most  rapidly  varying  solution,  by  using  the  standard 
predictor-corrector,  computer  time  may  be  wasted  in  integrating  the 
slower  varying  solutions  thru  use  of  the  short  step-size  h.  Thus, 
the  need  arises  for  a  generalized  predictor -corrector  algorithm  in 
which  the  step-sizes  may  be  different  in  each  eqi^tion.  To  this  end, 
we  wish  to  use  an  increment  h^  in  the  ith  equation  where  ^  »  m^h., 

i  =  2,  3,  ...,  N  and  m^  is  a  positive  integer. 

One  would  like  to  progress  to  the  point  +  h^  by  means  of  one 
predictor-corrector  step  in  the  first  equation,  m2  steps  in  the  sec¬ 
ond  equation,  m2  m^  steps  in  the  third  equation  ,  . . ,  m2  m^  . . .  m^. 

steps  in  the  Nth  equation.  That  is,  we  don't  wish  to  calculate  each 
fi  ■  y!  at  any  points  between  x^  +  Mh^  and  x^  +  (M+l)hj^,  i  »  1,  2, 

. . . ,  N. 


Suppose  for  the  moment  we  have  a  method  of  calculating 


y. (r.)  for  r.h. 
1  1  11 


h. 

j  *  1»  2,  ...,  T— ,  i  ■  1,  2,  ...,  N-1. 


Let  this  method  be  denoted  by  (*).  We  call  (*)  the  generalized  pre 
dictor  as  will  be  borne  out  in  the  following  discussion. 
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To  progress  from  to  we  use  the  following  procedure: 


N 

put  M  ■  1,  r^  *  IP  ,  i  «  1,  2,  N-1 

"i 

1.  predict  y^(M)  using  standard  predictor 

2.  predict  using  (*),  i  *  1,  2,  N-1 

3.  compute 

4.  correct  yjj(M)  using  standard  corrector 

5.  compute  using  corrected  y^jCM) 

put  M  =  M+1,  r^  =  2r^,  i  «^1,  2,  N-1 

times 

6.  compute  fj^_j(l)  using 

M  XnC^n) 

(b)  predicted  ^(1)  from  standard  predictor 
Cc)  y^(nij^  r^)  from  (*),  i  *  1,  2,  N-2 

7.  compute  fj^  ^(1)  using  corrected  yj^  ^(1) 

The  same  scheme  is  repeated  to  advance  from  x  h.,  ,  to  x  *  2h.,  , 
(using  (*)  advanced  one  step  for  Continuing  in  this  fashion 

then  finds  the  solution  at  x  +  m.,  ,  h  „  «  x  ♦  h..  We  are  now 

able  to  compute  a  corrected  value  for  y^  2  point.  When  we 

h. 

finally  reach  x^  hj,  we  have  evaluated  f^  at  only  jp  points. 


13 


i  ■  1,  2,  N.  For  a  simple  example  illustrating  the  above  pro¬ 
cedure  see  Appendix  A. 

Now  it  only  remains  to  determine  the  generalized  predictor  (*) 
for  calculating  the  predicted  value  of  where  y  is  the  inter¬ 

mediate  value  of  M,  0  <  Y  <  i*  To  do  this,  write  (20),  the  standard 
predictor,  in  the  form. 


PfCY)  '  aj(Y)  y^(M)  ♦  bj(Y)  y^(M-l)  ♦  c^(y)  y^(M-2) 


♦  dj(Y)  y.(M-3)  +  h.[ej(Y)  y!(M)  +  fj(Y)  y!(M-l)  (27) 

♦  gl(Y)  y|(M-2)  +  Kj(y)  y!(M-3)] 


where  the  constant  coefficients  are  to  be  determined  and  the  follow¬ 
ing  notation  is  used: 


Pi(Y)  *  y^(x^  +  Yb^) 
/^(M)  =  y^(x^  ♦  Mh.) 
y^(M)  =  y^(x^  Mh^) 


i  =  1 ,  2,  . . . ,  N 

i  =  1,  2,  ....  N  (28) 

% 

i  -  1 ,  2 ,  . . . ,  N 


By  the  same  technique  employed  before,  the  solution  set  is  found 

to  be 

aj(Y)  =  -  dj(Y)  -3ej(Y)  ♦  3Kj(y)  +  1/4  (y**  ♦  6y^  +  13y^  +  12  +  r) 
bj(Y)  =  -9dj(Y)  +  24Kj(y)  +  (y'^  +  4y^  +  4y^) 

Cj(y)  *  9dj(Y)  +  3e^(Y)  -  27Kj(y)  -  1/4  (Sy'*  +  22y^  ♦  29y^  +  12y) 
dj(Y)  =  dj(Y)  (29) 
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+  4y) 


CjCy)  =  CjCy) 

F^Cy)  =  -6d^(Y)  ♦  4ej(Y)  -  17K^(y)  -  (y^  +  5y^  ♦  8y^ 
gjCY)  “  6d^(Y)  ♦  ej(Y)  -  14Kj(y)  -  1/2  (y"*  +  4y^  +  Sy^  +  2y) 

The  corresponding  truncation  error  in  the  generalized  predictor 
is  given  by: 

E  =  12dj(Y)  -  4ej(Y)  ♦  40Kj(y) 

5  4  3  2 

*  (Y^  +  6y^  +  ISy"^  +  12y  ♦  4y)  -  (30) 


As  before,  the  standard  technique  for  improving  P^(y)  is  to  add 


E(y) 


‘■i) 


t31) 


F:(y)»  where  E,  is  E  at  y  *  1.  E-  is  from  (23),  P.  is  the  pre- 
X  1  Y  C  A  X 

dieted  value  of  the  previous  step  and  C  is  the  corrected  value  of 
the  previous  step. 

I 

Similarly,  the  truncation  error  in  standard  corrector  is  reflected 
in 

Et  (M.1)  =  (P^(M.l)  -  C^(M.l)) 

c  c  p  (32) 


where  E^  and  E^  are  the  respective  truncation  errors  in  the  standard 

predictor-corrector.  By  calculating  E.j.  (M+1),  the  validity  of  the 

c 

answer  found  by  using  the  standard  predictor-corrector  may  be  checked. 
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However,  the  task  of  accounting  for  the  truncation  error  in  the 
steps  using  the  generalized  predictor  is  somewhat  different.  Since 
only  predicted  values  from  the  generalized  predictor  are  available 
at  intermediate  steps,  the  truncation  error  in  the  generalized  pre¬ 
dictor  is  reflected  in 


Et  (y)  •  f ^  i  »  1.  2,  ...,  N-l  (33) 

^c  "  ^1  "  ^ 


where  is  at  y  *  i*  and  are  the  predicted  and  corrected 

values,  respectively,  obtained  at  the  current  step.  This  value  must 
be  used  throughout  the  current  stage. 


Analysis  of  Predictor-Corrector  Methods 


One  basic  shortcoming  in  all  predictor -corrector  methods  is  that 
they  require  some  other  procedure  to  obtain  enough  values  so  as  the 
method  can  be  started.  In  the  case  of  the  generalized  predictor- 
corrector,  some  procedure  must  be  used  to  generate  the  first  3(m^) 

steps  in  the  ith  equation.  One  method  of  starting  the  proceudre  is 
the  Runge-Kutta  method. 

One  distinct  disadvantage  of  the  predictor -corrector  methods 
is  the  difficulty  in  halving  the  step-size  h.  By  utilizing  the  gen¬ 
eralized  predictor-corrector  some  of  these  difficulties  can  be  sur¬ 
mounted.  Still,  sometimes  halving  one  of  the  increments  is  required, 
necessitating  interpolation  or  even  restarting  the  current  step. 

Stability  of  the  predictor-corrector  algorithms  was  not  discus¬ 
sed  here  due  to  the  lengths  that  one  is  required  to  go  to  an  ade¬ 
quate  discussion.  Suffice  it  to  say,  the  system  of  differential 
equations  to  be  solved  must  be  carefully  examined  to  determine  what 
odd  characteristics  exist  and  then  the  proper  type  of  predictor-cor¬ 
rector  algorithm  chosen  to  give  the  desired  stability  and  accuracy. 
For  a  further  discussion  see  [10] . 

Predictor-corrector  methods  in  general^  are  much  faster  than  the 
Runge-Kutta  methods  and  in  particular,  use  of  the  generalized  predic¬ 
tor  will  reduce  the  function  evaluations  per  step  and  hence  further 
reduce  computer  time  used. 
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Block  methods  can  also  be  applied  to  predictor>corrector  schemes 
to  improve  their  accuracy  and  efficiency.  For  a  more  complete  dis¬ 
cussion  of  block  methods,  see  [17]. 


Hybrid  Method 


Hybrid  methods  differ  from  those  previously  described  in  that 
in  addition  to  using  previous  information  from  the  M,  (M-l)st,  (M-2)nd, 
etc.  steps,  an  outside  method  is  used  to  calculate  information  at 
some  (M-Y)th  step,  0  <  y  <  1,  and  this  information  is  also  incorpoj  .»tsd 
into  the  multistep  procedure. 

In  certain  cases  the  added  information  tends  to  stabilize  the 
formula  in  which  it  is  used  and  makes  halving  the  step-size  an  easier 
task. 


However,  as  pointed  out  in  Gear  [7],  situations  arise  in  which 
the  use  of  the  added  information  introduces  much  more  error  into  the 
solution  of  the  system  than  is  desirable.  For  a  more  thorough  dis¬ 
cussion  of  hybrid  methods,  see  [7]  and  [9] . 


SUMMARY 


The  Runge-Kutta  algorithm  can  be  used  for  the  solution  of  most 
systems  of  differential  equations.  Since  only  information  from  one 
previous  step  is  required,  the  Runge-Kutta  methods  have  favorable 
stability  characteristics  and  the  process  of  altering  the  step-size 
is  an  easy  task.  Consequently,  Runge-Kutta  methods  are  frequently 
used  as  a  starting  procedure  for  other  methods  such  as  multistep  methods. 
However,  the  number  of  function  evaluations  per  step  required  by  the 
Runge-Kutta  methods  requires'utilizing  more  computer  time  than  multi- 
step  methods  and  since  no  measure  of  the  truncation  error  is  avail¬ 
able  during  the  calculation  procedure,  approximate  or  comparison  methods 
are  required  in  the  program  to  reflect  this  error. 

Multistep  methods  require  fewer  function  evaluations  per  step 
and  hence,  use  less  computer  time.  The  truncation  error  can  be  meas¬ 
ured  effectively  during  the  calculation  procedure,  thereby  exhibit¬ 
ing  the  error  more  correctly  than  do  Runge-Kutta  methods.  The  sta¬ 
bility  of  the  predictor-corrector  method  should  be  investigated  thor¬ 
oughly  to  determine  what  regions  of  integration  are  required  for 
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efficient  use  of  these  methods.  Although  a  smaller  step«size  is  gen¬ 
erally  required  for  predictor-corrector  methodsras  compared  to  the  Runge- 
Kutta  method,  the  use  of  the  generalized  predictor-con ector  method 
may  alleviate  many  problems  found  in  multistep  methods  and  may  use 
the  least  computer  time  of  all  the  methods  described. 

The  fairly  attractive  methods  briefly  described  in  the  final 
section  point  out  other  method^ that  stem  from  the  Runge-Kutta  and 
the  multistep  methods,  that  may  be  investigated  for  use  in  a  ballistic 
rocket  trajectory  program  due  to  their  composite  of  the  desirable 
characteristics  of  the  Runge-Kuttr  and  multistep  methods. 

At  the  present  time  the  Atmospheric  Sciences  Office  at  White 
Sands  Missile  Range  is  programming  the  generalized  predictor -correc¬ 
tor  algorithm  for  use  in  the  six-degree  of  freedom  ballistic  rocket 
model.  It  is  felt  that  this  algorithm  can  best  accomplish  the  desired 
task  with  a  minimum  of  computer  time  and  with  accuracy  comparable 
to  that  of  Runge-Kutta  methods  used  currently. 
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APPENDIX  A 


An  Example  Using  the  Generalized  Predictor-Corrector  Method 

Suppose  we  have  the  following  system  of  first  order  differential 
equations 

Xl'  -  (X.  72 .  Xj) 

X2'  =  ^2  ^1*  ^2’  ^3^ 

^3'  *  ^3  ^1*  ^2*  ^3^ 


where  72  varies  three  times  as  rapidly  as  y^  with  respect  to  x,  and 
yj  varies  twice  as  rapidly  as  72  with  respect  to  x.  We  would  like 
to  find  the  solution  of  (34)  at  the  point  x^  +  hj.  Therefore,  let 
hj  =  3h2  and  h2  =  2hj  where  h^  is  the  increment  in  the  ith  equation. 
This  means  that  :a^  =  3  and  m^  =  2.  The  following  figure  will  illus¬ 
trate  this  situation. 


In  finding  the  solution  of  (34)  at  x^  +  h^  we  will  calculate 

at  x^  +  hj  only,  £2  at  x^  +  h2  and  x^  +  2h2  only  and  at  x^  +  Kh^, 

K  =  1,  2,  ...,  6,  and  use  the  generalized  predictor  for  finding  the 
y's  between  these  points. 
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The  detailed  procedure  is  the  following: 

1.  put  r^^  a  ^  ,  i  »  1,  2 

"i 

2.  a.  predict  (1)  using  the  standard  predictor 

b.  predict  y^  (r^)  =  y^  (1/6)  using  (*) 

c.  predict  (r2)  =  y2  (1/2)  using  (*) 

3.  compute  (1),  correct  y^  (1)  using  the  standard  corrector 

and  then  compute  a  final  value  of  f^  (1) 

4.  put  r.  a  2r.  =  m-  r.,  i  =  1,  2 

5.  a.  predict  y^  (2)  a  y^  (m^)  using  the  standard  predictor 

b.  predict  y^  (r^)  =  y^  (2/6)  using  (•) 

c.  predict  y2  (r2)  =  y2  (1)  using  the  standard  predictor 

6.  compute  f^  (2),  correct  y^  (2)  using  the  standard  corrector, 
and  then  compute  a  final  value  of  f^  (2) 

7.  compute  f2  (1)  using  5,  correct  y^  (1)  using  the  standard 
corrector,  and  then  compute  a  final  value  of  f2  (1) 

In  the  second  equation  we  are  now  at  the  point  x  +  h2.  We  now  ad¬ 
vance  (*)  one  step  when  using  it  in  the  second  equation  and  it  will 
be  understood  that  when  we  say  "predict  y2  (r2),  0<r2<l,  using  (*)" 

h’e  mean  we  are  predicting  the  value  of  y2  between  +  h2  and  x  +  2h 

hj 

8.  put  r^  =  3r^  =  (m^  +  1)  rj^,  r2  =  jp 

9.  a.  predict  y^  (3)  using  the  standard  predictor 

b.  predict  y^  (r^,)  =  y^  (3/6)  using  (*) 

c.  predict  y^  (r^)  »  y^  (1/2)  using  (*) 
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10.  c(»pute  (3),  correct  (3)  using  the  standard  corrector 
and  then  compute  a  final  value  of  f^  (3) 

11.  put  r^  «  4rj  =  2inj  r^,  T2  *  2t2  •  m,  T2 

12.  a.  predict  y^  (4)  ■  y^  (2inj)  using  the  standard  predictor 

b.  predict  y^  (r^)  »  y^  (4/6)  using  (*) 

c.  predict  y2  (r2)  *  y2  (1)  using  the  standard  predictor 

13.  COTipute  fj  (4),  correct  y^  (4)  using  the  standard  corrector, 
and  then  compute  a  final  value  of  f^  (4) 

14.  compute  f2  (1)  using  12,  correct  y2  (1)  using  the  standard 
corrector  and  then  compute  a  final  value  of  f2  (1) 

Again  we  advance  (*)  one  step  when  using  it  in  the  second  equation. 

We  will  now  compute  y2  (r2),  o<r2<l  between  2h2  and  ♦  3h2. 

^3° 

15.  put  rj  =  5rj  =  (2mj  *  1)  r2  »  jj- 

16.  a.  predict  y^  (5)  using  the  standard  predictor 

b.  predict  y^  (r^)  *  y^  (5/6)  using  (*) 

c.  predict  y2  (r2)  «  y2  (1/2)  using  (*) 

17.  compute  fj  (5),  correct  y^  (5)  using  the  standard  corrector, 
and  then  compute  a  final  value  of  f^  (5) 

18.  put  rj  *  6rj  =  m^  m2  r^,  r2  ®  2r2  =  m^  r2 

19.  a.  predict  y^  (6)  =  y^  (m^  m2)  using  the  standard  predictor 

b.  predict  y^^  (r^)  =  y^  (1)  using  the  standard  predictor 

c.  predict  y2  (r2)  =  y2  (1)  using  the  standard  predictor 

20.  from  19,  compute  f^  (1),  correct  y^  (1)  using  the  standard 
corrector,  and  then  compute  a  final  value  of  f^^  (1),  i  =  1,2,3. 
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We  have  reached  our  solution  at  h^  by  means  of  one  standard  pre¬ 
dictor-corrector  step  in  the  first  equation,  ®2  ®  ^  standard  predictor 
corrector  steps  in  the  second  equation,  and  =  6  standard  pre¬ 

dictor-corrector  steps  in  the  third  equation. 
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